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Surface reflectance adjusted for atmospheric effects is a primary input for land cover change detection 
and for developing many higher level surface geophysical parameters. With the development of 
automated atmospheric correction algorithms, it is now feasible to produce large quantities of surface 
reflectance products using Landsat images. Validation of these products requires in situ measurements, 
which either do not exist or are difficult to obtain for most Landsat images. The surface reflectance 
products derived using data acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS), 
however, have been validated more comprehensively. Because the MODIS on the Terra platform and 
the Landsat 7 are only half an hour apart following the same orbit, and each of the 6 Landsat spectral 
bands overlaps with a MODIS band, good agreements between MODIS and Landsat surface reflectance 
values can be considered indicators of the reliability of the Landsat products, while disagreements may 
suggest potential quality problems that need to be further investigated. Here we develop a system 
called Landsat-MODIS Consistency Checking System (LMCCS). This system automatically matches 
Landsat data with MODIS observations acquired on the same date over the same locations and uses 
them to calculate a set of agreement metrics. To maximize its portability, Java and open-source libraries 
were used in developing this system, and object-oriented programming (OOP) principles were followed 
to make it more flexible for future expansion. As a highly automated system designed to run as a stand- 
alone package or as a component of other Landsat data processing systems, this system can be used to 
assess the quality of essentially every Landsat surface reflectance image where spatially and temporally 
matching MODIS data are available. The effectiveness of this system was demonstrated using it to 
assess preliminary surface reflectance products derived using the Global Land Survey (GLS) Landsat 
images for the 2000 epoch. As surface reflectance likely will be a standard product for future Landsat 
missions, the approach developed in this study can be adapted as an operational quality assessment 
system for those missions. 

© 2011 Elsevier Ltd. All rights reserved. 


1. Introduction 

Surface reflectance is the most basic remotely sensed earth 
surface parameter in the solar reflective wavelengths (Vermote 
et al., 2006). It provides a primary input for retrieving essentially 
most higher level surface geophysical parameters, including vege- 
tation indices (Huete et al, 2002), albedo (Fang et al, 2007), leaf 
area index (Myneni et al., 1997), and photosynthetically active 
radiation (PAR) (Liang et al., 2006, 2007). While both raw digital 
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numbers (DN) and calibrated top-of-atmospheric (TOA) reflec- 
tance have been widely used in land cover mapping and change 
detection, use of surface reflectance can improve the results from 
such analysis (Song et al., 2001). 

With the launch of the first Landsat in 1972, a series of Landsat 
satellites have produced large quantities of images useful for land 
cover and change studies and other earth science applications 
(Townshend et al M 1991; Goward and Williams, 1997; Goward 
et al., 2006; Huang et al., 2009a, b). Despite this near 40-year 
history of the Landsat missions, however, standard Landsat 
imagery products have been provided in DN or higher level 
products (e.g., TOA (Roy et al., 2010)), not surface reflectance. 
This is due in part to many challenges to performing atmospheric 
correction operationally on Landsat images. Recently, the 
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atmospheric correction algorithm developed for the Moderate 
Resolution Spectroradiometer (MODIS) (Vermote et al., 2002) was 
adopted for use with Landsat data, and was implemented as part 
of the Landsat Ecosystem Disturbance Adaptive Processing 
System (LEDAPS) developed at NASA Goddard Space Flight Center 
(GSFC) (Masek et al., 2006). As an adaptation of the MODIS 
Adaptive Processing System (MODAPS) (Justice et al, 2002) for 
processing Landsat data, the LEDAPS allows rapid processing of 
large quantities of Landsat images to produce surface reflectance 
products from the raw radiometry. It has been used to produce a 
surface reflectance record consisting of over 2000 Landsat images 
over North America (Masek et al., 2006). After expanding its 
computing and storage capacity, this system is being adopted for 
creating global land surface reflectance products using Landsat 
images assembled through the joint National Aeronautics and 
Space Administration (NASA)-U.S. Geological Survey (USGS) 
Global Land Survey (GLS) program, which provide wall-to-wall, 
orthorectified Landsat coverage of Earth’s land area for four 
epochs centered around 1975, 1990, 2000, and 2005 (Gutman 
et al., 2008). The LEDAPS atmospheric correction algorithm has 
also been adopted by other research groups that have needs for 
processing large quantities of Landsat images, and is being 
considered a candidate algorithm for producing surface reflec- 
tance products for the next Landsat system, the Landsat Data 
Continuity Mission (LDCM) (Irons and Masek, 2006). 

The usefulness of a data product depends to a large degree on 
its quality and accuracy. Knowledge of the uncertainties of a 
product is required in order to understand the impact of such 
uncertainties on downstream applications. While some surface 
reflectance products derived from Landsat Thematic Mapper (TM) 
and Enhanced Thematic Mapper Plus (ETM + ) images have been 
validated using ground observations (Vermote et al., 1997a; Liang 
et al., 2002; Masek et al., 2006; Ouaidrari and Vermote, 1999), the 
applicability of this ground-based validation approach is 
constrained by the availability of ground observations, including 
measurements from AERONET sites (Wang et al., 2009). For large 
area applications using Landsat images, the required images are 
typically acquired in dates that most likely differ by months, 
seasons, or even years. As a result, even when ground observa- 
tions are available for validating the reflectance value of some of 
those images, the conclusions reached may not be applicable to 
other images. In addition to comparison with ground observa- 
tions, the quality of a product can also be assessed with respect to 
their intended performance by comparison with other existing 
products (Roy et al., 2002). Since MODIS surface reflectance 
products have been validated comprehensively (Vermote et al., 
2002; Liang et al., 2002; Fang et al., 2004; Kotchenova and 
Vermote, 2007; Vermote and Kotchenova, 2008a), they may be 
used as a reference to assess the quality of Landsat surface 
reflectance products. The MODIS on the Terra platform follows 
the same orbit as the Landsat 7. It crosses the equator at about 
10:30 AM local solar time, roughly 30 min later than the Landsat 
7 (Vermote and Kotchenova, 2008b). For each of the 6 Landsat 
multispectral bands, there is a MODIS band with a spectral range 
well within that of the corresponding Landsat band (Landsat 
7 Science Data Users Handbook, 2009) (see Table 1). Comparisons 
between MODIS and Landsat surface reflectance products were 
performed for limited sites in previous studies, which revealed 
that the differences in bandwidth between the two instruments 
only resulted in tolerable differences in reflectance values 
(Vermote et al., 2002; Masek et al., 2006). The moderate resolu- 
tion MODIS data are coarser than the 30 m resolution Landsat 
data, but they provide global, near-daily repeat coverage, 
compared to the 16-day repeat cycle of Landsat 7. In addition, 
MODIS Terra data products have been made available globally 
since 2000; they allow quality assessment (QA) of essentially all 


Table 1 

The 6 multispectral bands of the Landsat 7 and the corresponding MODIS bands. 


Landsat ETM + 
band 

ETM+ bandwidth 
(nm) 

MODIS 

band 

MODIS bandwidth 
(nm) 

1 

450-520 

3 

459-479 

2 

530-610 

4 

545-565 

3 

630-690 

1 

620-670 

4 

780-900 

2 

841-876 

5 

1550-1750 

6 

1628-1652 

7 

2090-2350 

7 

2105-2155 


surface reflectance products that could be produced using Landsat 
7 images acquired since 2000. 

The purpose of this paper is to develop a Landsat-MODIS 
Consistency Checking System (LMCCS) for systematically check- 
ing the consistency between Landsat 7 surface reflectance pro- 
ducts and MODIS Terra products. High levels of agreements 
between the two sets of products are supporting evidence that 
the surface reflectance values in the Landsat products are com- 
parable with those in the MODIS data, while obvious discrepan- 
cies between them likely indicate data quality problems with one 
or both sets of products that warrant further examination. In this 
approach, object-oriented programming (OOP) principles have 
been implemented in the system design to support further model 
extension. The system development has been facilitated using 
open-sources libraries. This system could be used as a stand-alone 
system or together with a Landsat processing system like the 
LEDAP. This paper describes the major considerations in design- 
ing and implementing this system, and highlights some of the 
issues revealed by this system in evaluating preliminary surface 
reflectance products derived using the GLS 2000 Landsat data set. 

2. Overview of MODIS and Landsat surface reflectance 
algorithms 

2.1. MODIS surface reflectance algorithm 

The MODIS atmospheric correction scheme is based on the 
Second Simulation of a Satellite Signal in the Solar Spectrum 
Vector Code (6S) radiative transfer code (Vermote et al., 1997b). It 
calculates surface reflectance by accounting for the effect of 
gaseous absorption, molecules, and aerosol scattering. Some 
limitations of the algorithm are the use of a Lambertian model 
especially at shorter wavelength (e.g., 0.55 pm) (Wang et al., 
2010) and the use of the dark target approach for aerosol retrieval 
that does not produce results over bright targets like deserts (Levy 
et al., 2007). Detailed descriptions of the MODIS surface reflec- 
tance algorithm can be found in many previous publications 
(Vermote et al., 1997a, 2002; Vermote and Kotchenova, 2008a). 

Several MODIS surface reflectance products produced at differ- 
ent spatial resolutions using different temporal aggregation or 
compositing schemes are available (Justice et al., 2002). The MODIS 
Terra daily gridded surface reflectance product (MOD09GA) avail- 
able for download from NASA’s Warehouse Inventory Search Tool 
(WIST, https://wist.echo.nasa.gov/) is used in this study. This 
product has the Sinusoidal projection and provides reflectance 
value at the 500 m spatial resolution for MODIS bands 1-7, as well 
as observation and geolocation statistics. 

2.2. Landsat surface reflectance algorithm 

The Landsat surface reflectance algorithm is described in 
Masek et al. (2006) and is based on the MODIS approach 
(Vermote et al., 2002). It uses the 6S radiative transfer code to 
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compute the transmission, intrinsic reflectance, and spherical 
albedo for relevant atmospheric constituents (Vermote et al, 
1997b), and to calculate surface reflectance by compensating for 
atmospheric scattering and absorption effects on the TOA reflec- 
tance. The relevant atmospheric constituents include gases, 
ozone, water vapor, and aerosols. Ozone concentration was 
derived from the Total Ozone Mapping Spectrometer (TOMS) 
aboard the Nimbus-7, Meteor-3, and Earth Probe platforms as 
well as from NOAA’s Television Infrared Observation Satellite 
Program (TIROS) Operational Vertical Sounder (TOVS) ozone data 
when TOMS data were not available. Column water vapor was 
taken from NOAA National Centers for Environmental Prediction 
(NCEP) reanalysis data (available at http://dss.ucar.edu/datasets/ 
ds090.0/). Digital topography (1 km GTopo30) and NCEP sea level 
pressure data were used to adjust Rayleigh scattering to local 
conditions. Aerosol optical thickness (AOT) was directly derived 
from the Landsat image using the dark, dense vegetation method 
of Kaufman et al. (1997). 


3. Considerations for comparing Landsat and MODIS surface 
reflectance products 

3.1. Agreement measures 

Since the Landsat 7 and the Terra MODIS follow the same orbit, 
their equator crossing times differ by only about 30 min (i.e., 
equatorial crossing times from 10:00 to 10:15 AM for Landsat 
7 and 10:30 for Terra), and each Landsat multispectral band has a 
MODIS band with a similar bandwidth (Landsat 7 Science Data 
Users Handbook, 2009) (Table 1), the Landsat reflectance values, 
when aggregated spatially using the MODIS pixel coverage, 
should be similar to the MODIS reflectance values for the same 
location. For a Landsat image acquired at date (t), its reflectance 
values can be aggregated using the MODIS pixel coverage accord- 
ing to 

C t =5>jxp[), (1) 

where C t is an aggregated reflectance value, i is spatial index 
(location) of Landsat pixels within the concerned MODIS pixel 
coverage, l\ is surface reflectance of a Landsat pixel, and p\ is area 
percentage of that pixel in the MODIS pixel coverage. Since 
Landsat and MODIS surface reflectance values are expected to 
be very close, their relationships can be modeled using a simple 
linear regression as follows to examine the agreement between 
Landsat (C t ) and MODIS surface reflectance ( M(x if y Jf s,t )): 

M(x it y jt s,t ) = a x C(Xi,y jt s,t ) + b, (2) 

here (x i} yj) (ij e [1 ,n]) is a given pixel location for both Landsat and 
MODIS images. Suppose we have n pairs input of C(x it yj,s,t ) and 
M(Xi,yj,s,t) and each pair is acquired on the same date, s and t 
represent the spectra band and acquisition date for both MODIS 
and ETM+ data, a represents the linear slope, and b represents 
the offset of linear relationship between MODIS and ETM + 
surface reflectance, a and b can be calculated to represent the 
differences between the surface reflectance pairs using linear 
regression, and represent the line that has minimum error value 
from all input data regardless of the correlation and difference 
between the two input values. In order to evaluate the correlation 
and differences between MODIS and Landsat data, R 2 (see Eq. (3)) 
and root mean square deviation (RMSD) (see Eq. (4)) are intro- 
duced into the approach 

r2 1 E i j(M(x,-,y j ,t)-C(x i .y J ,t)) 2 (3) 




3.2. Spatial considerations 

Use of the above agreement measures as indicators of the 
consistency between Landsat and MODIS reflectance values 
requires that the ground footprints of the two sets of values 
match. Matching the footprints of Landsat and MODIS observa- 
tions, however, can be complicated by several issues, including 
data coordinate reference system (CRS), residual misregistration 
errors (Rojas et al., 2002), sensor’s point spread function (PSF) 
(Vermote et al., 1997a; Huang et al., 2002), adjacency effect 
(Ouaidrari and Vermote, 1999), and view angle effect (Schaaf 
et al., 2002). These issues need to be addressed properly in order 
to minimize their impact on the agreement measures calculated 
using Eqs. (l)-(4). 

3.2.1. Reprojection of MODIS data 

The MODIS data are provided in the Sinusoidal projection, but 
the Landsat images have the Universal Transverse Mercator 
(UTM) projection. Because the spatial resolution of the MOD09GA 
data is 500 m, roughly 16.6 times that of the Landsat data, the 
reflectance values in the Landsat data need to be aggregated 
spatially to match the MODIS pixels. In order to match the 
coverage of each MODIS pixel with correspondent Landsat pixels 
more accurately, the four corners of each MODIS pixel can be used 
to draw a geospatial polygon to represent the boundary of that 
MODIS pixel, which can then be reprojected using the projection 
of the Landsat data (to be discussed in detail in Section 4.2). Only 
Landsat pixels located within the boundary of that MODIS pixel 
should be used to calculate an average reflectance value, which 
will be compared to the value of that MODIS pixel (see Fig. 1). 

3.2.2. Use of homogeneous regions 

The Landsat-MODIS relationship as described in Eq. (2) likely 
will degrade over nonhomogeneous areas, because the impact of 
residual misregistration errors (Rojas et al., 2002), adjacency 
effect (Ouaidrari and Vermote, 1999), and sensor’s PSF is more 
significant in such areas (Vermote et al., 1997a; Huang et al., 
2002). The residual geolocation errors in the GLS Landsat data 
could be up to 50 m (lcr) (Tucker et al., 2004), and the 500 m 
resolution MOD09GA data geolocation has also been improved to 
50m (lcr) at nadir (Wolfe et al., 2002). While in Eq. (1) the 
Landsat pixels are aggregated using a simple spatial average 
method, the actual MODIS PSF has higher weighting factors at 
the center of a MODIS pixel and lower weights toward pixel 
edges, and has non-zero weights for immediately neighboring 
areas outside the footprint of that pixel (Rojas et al., 2002). In 
nonhomogeneous areas considerable differences may exist 
between the values derived using the MODIS PSF and a simple 
spatial average filter (Huang et al., 2002). The impact of the 
residual misregistration errors and the sensor’s PSF can be 
minimized using samples selected from homogeneous areas. 

Since the MOD09GA data have a geolocation accuracy of less 
than one MODIS pixel (Wolfe et al., 2002), we use a window with 
size of 3 x 3 MODIS pixels to identify homogeneous regions at the 
MODIS resolution. Specifically, for each pixel a value range is 
calculated as the difference between the maximum and the 
minimum values of the 9 MODIS pixels in the 3x3 window 
centered at that MODIS pixel. A MODIS pixel is considered in a 
homogeneous region if the value range is less than a predefined 
threshold value. We developed the value range threshold based 
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Fig. 1. A polygon defined by the four corners of a MODIS pixel delineates the boundary of that pixel. The polygon is reprojected using the projection used by a Landsat 
image (i.e., UTM) before it is overlaid on the Landsat image to determine the Landsat pixels that are located within the boundary of that MODIS pixel. 


Table 2 

Value range threshold values used to determine homogeneous regions at both 
Landsat and MODIS resolutions. 


MODIS band 

Landsat band 

Threshold (absolute reflectance value) 

3 

1 

0.03 

4 

2 

0.03 

1 

3 

0.03 

2 

4 

0.06 

6 

5 

0.03 

7 

7 

0.03 


on our observations at homogeneous Landsat pixels for the 
6 Landsat bands (see Table 2). Similarly, the value range of all 
Landsat pixels located within a MODIS pixel coverage is used to 
determine whether the MODIS pixel is also homogeneous at the 
Landsat resolution. Only samples that are homogeneous at both 
the Landsat and the MODIS resolutions are selected. These 
samples are then divided into 10 equal percentiles using a 
cumulative histogram, and 20% of the samples in each percentile 
are used to calculate the agreement measures described in 
Section 3.1. 

3.2.3. View angle control 

In mid- to high-latitude regions, MODIS can produce more 
than one observation per day in certain areas (Wolfe et al., 2002). 
In those areas the pixels in the MOD09GA product do not 
necessarily have the same view zenith angle as the Landsat image 
acquired on the same day. To avoid the potential impact of 
different view zenith angles on the reflectance values, only MODIS 
pixels having view zenith angles within the view zenith angle 
range of the Landsat (i.e., + 7.5° from nadir) are used in the 
comparison. 

3.3. Temporal considerations 

Because the MODIS on the Terra platform is only about 30 min 
later than the Landsat 7, the approach assumes that the atmo- 
spheric condition does not change much in 30 min; the two 


Table 3 

The three bits used to indicate cloud state and shadow presence in the MODIS 
1 km resolution QA band (Vermote and Vermeulen, 1999). 


Bit no. 

Parameter name 

Bit combination 

State 

0-1 

Cloud state 

00 

Clear 



01 

Cloudy 



10 

Mixed 



11 

Not set 

2 

Cloud shadow 

1 

Yes 



0 

No 


instruments should see roughly the same atmospheric condition 
when they pass the same location each day (Vermote et al., 2002). 
When clouds are present in their view, however, the two instru- 
ments may find those clouds and their shadows at different 
locations, because clouds can move a lot in 30 min. To minimize 
the impact due to cloud movement on the comparison results, we 
exclude both cloud and shadow pixels from the comparison. In 
the MODIS data, cloud and shadow are identified using the 
internal cloud and shadow masks from the MODIS QA band, 
which is a 16 bits band with 1 km spatial resolution. Two of the 
16 bits are used to represent cloud state; one bit indicates cloud 
shadow status (Vermote and Kotchenova, 2008b) (see Table 3). 
Since minimal cloud cover was one of the considerations in 
selecting the GLS images (Tucker et al., 2004), the LMCCS tool 
does not check the cloud cover in the Landsat images. 

3.4. Key processing steps 

Based on the above considerations, four major steps are 
needed in order to perform a Landsat-MODIS comparison for 
each spectral band: 

1. Initial sample selection. For each Landsat image and its sim- 
ultaneously acquired MODIS image, the comparison is per- 
formed based on samples. Each sample is a surface reflectance 
value pair collected from both Landsat and MODIS images at 
the same coordinate. The MODIS value is retrieved from the 
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MODIS pixel at the coordinate, and the Landsat value is 
retrieved by aggregating the Landsat pixels within the MODIS 
pixel coverage (see next entry). The sample coordinates could 
be selected either randomly or with the use of a systematic 
sampling method within the overlaid region of the two 
images. In the system developed in this study, the sample 
coordinates are selected systematically by picking one out of 
3 MODIS pixels in both the horizon and vertical directions. 

2. Spatial aggregation of Landsat data. As discussed above, in order 
to match a MODIS pixel with the Landsat pixels, a polygon 
is generated to represent the MODIS pixel coverage under 
MODIS Sinusoidal projection, which is reprojected using the 
UTM projection used by the Landsat image. The Landsat 
pixels located within the reprojected polygon are used to 
calculate an average surface reflectance value that will be 
compared with the MODIS reflectance value at that sample 
coordinate. 

3. Sample filtering. The initial samples are filtered to exclude 
nonhomogeneous pixels, pixels contaminated by cloud and 
shadow, MODIS pixels with view zenith different than the 
Landsat, etc. The samples used in the final comparison are 
selected to represent the spectral range of the concerned images. 


4. Calculation of agreement measures. The samples retained after 
the filtering process in step 3 are used to calculate the 
agreement measures according to Eqs. (3) and (4). 


4. Implementation of the LMCCS 

4.1. System overview 

The LMCCS consists of two modules, i.e., Geospatial Image 
Metadata Administration (GIMA) and Surface Reflectance Consis- 
tency Checking (SRCC) module (Fig. 2). The GIMA module ingests 
metadata from collected MOD09GA images to facilitate quickly 
identifying simultaneously acquired MODIS image for the 
Landsat-MODIS comparison. The SRCC module implements the 
comparison discussed in Section 3. Open-source geospatial 
libraries, e.g., GeoTools (http://www.geotools.org), PR0J.4 (http:// 
trac.osgeo.org/proj/), and JFreeChart (http://www.jfree.org/jfree 
chart), are used in system development. To add flexibility, the 
system was designed to run as a stand-alone package or as a 
component of other Landsat data processing systems. 


MODIS Surface 
Reflectance Images 



MODIS Metadata 

(Spatial Coverage, 
Acquired date. . .) 


Landsat Metadata 

(Spatial Coverage, 
Acquired date. . .) 


Landsat -MODIS Consistency Checking System (LMCCS) 


Geospatial Image Metadata Administration (GIMA) 


Ingest MODIS images 
metadata 


Identify MODIS 
images 



Landsat Surface 
Reflectance Image 


Request 


Perform consistency 
checking 


Output charts & 
metrics 


Surface Reflectance 
Consistency Checking (SRCC) 


System Integration for Automated Environment 


Landsat Data 
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Module/System 



MODIS Image 
Metadata Database 


-> / 


/ r ') 

ft. 


Agreement 
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Fig. 2. A flowchart of the Landsat-MODIS Consistency Checking System (LMCCS). 
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4.2. Geospatial Image Metadata Administration (GIMA) 

Each GLS epoch includes about 9000 Landsat images, and 
roughly the same amount of MODIS images need to be collected 
to perform the Landsat-MODIS comparison. For each Landsat 
image, identifying its simultaneously acquired MOD09GA images 
requires spatial and temporary matching analysis within the 
available MOD09GA images (see Section 4.2.2), and it would be 
inefficient to perform the analysis without database support, 
especially when processing thousands of Landsat images. In 
addition, it would be very challenging to manually track the 
image acquisition process and manage such a large amount of 
images. To address this problem, we built a database to manage 
the acquired MODIS images and provide an efficient way for 
identifying MODIS images. 


4.2.1. Geospatial powered MODIS Image Metadata Database 

The metadata for each MOD09GA image includes nonspatial 
components (e.g., acquisition date, image storage path, etc.), and 
spatial components (e.g., geospatial coverage). In order to support the 
storage and operations of geospatial metadata items, the MODIS 
Image Metadata Database is built based on the PostGIS package to 
leverage its geospatial data storage and operation capabilities. PostGIS 
extends the PostgreSQL database and supports geospatial data types 
and functions defined in OGC Simple Features Specification for SQL 
(Egenhofer, 1990). In addition to the columns for those nonspatial 
metadata items, a geospatial polygon-type column is added to store 
the spatial coverage of each image and Generalized Search Tree (GiST) 
geospatial indexing is created on the column to enhance performance 
of geospatial query. Unfortunately, the Sinusoidal projection used by 
the MODIS data is not supported by the current version of PostGIS. To 
address this problem, we added this projection to PostGIS using 
PROJ.4 and Well-Known Text (WKT) definition syntax: 

PROJ.4: +proj = sinu + lon_0 = 0 +x_0 = 0 +y_0 = 0+a= 6371 007. 
181 + b =6371007.1 81 + units=m 

WI<T: PROJCS[”MODIS_Sinusoidal”, GEOGCS[”GCS_Sphere”,[”D_- 
Sphere”,SPHEROID[”Sphere”, 6371007.181, 0.0]], 
PRIMEM[”Greenwich”, 0. 0 ], UNIT] ’De- 
gree ”,0.01 74532925 1 99433 ] ], 

PROJECTION[”Sinusoidal”],PARAMETER[”False_Easting”,0.0], 


PARAMETER] ”False_Northing”, 0. 0 ], PARAMETER] ”Cen tra l_ 

Meridian”, 0.0], UNIT[”Meter”, 1-0]] 

The GIMA allows users to retrieve metadata from the 
MOD09GA images in the HDF-EOS format (Wei et al., 2007), and 
ingest these retrieved metadata items into the MODIS Image 
Metadata Database. Because MODIS data are provided using a 
tiling system that does not change over time, and the spatial 
coverage of each tile predefined and known, the spatial coverage 
polygon of each MODIS image is derived based on its path/row 
number in the MODIS tiling system (see Fig. 3A). 

4.2.2. Identification of simultaneously acquired MODIS images 

Identifying simultaneously acquired MODIS images is compli- 
cated by handling the different projections and tiling systems 
of MODIS and Landsat data. The MOD09GA data are distributed 
by NASA WIST in MODIS Sinusoidal projection and MODIS tiling 
system (see Fig. 3 A). Landsat TM/ETM+ data are distributed 
in UTM projection and USGS Worldwide Reference System-2 
(WRS-2) tiling system (see Fig. 3B). 

For each Landsat image, the MODIS images that were acquired 
on the same date as the Landsat image and encompass or 
intersect with that image need to be identified. This can be 
achieved by querying the geospatial MODIS Image Metadata 
Database. One complication is that the MODIS images have the 
Sinusoidal projection while the Landsat images have the UTM 
projection. The spatial coverage of each Landsat image needs be 
converted to the MODIS Sinusoidal projection before it can be 
overlaid on the MODIS images. 

The mathematic construction of the MODIS Sinusoidal projec- 
tion at coordinates (x,y) is defined as 

x = R(2-2 0 )cos(p y = Rep, (5) 

where R is the radius of the Earth defined as a sphere 
(6,371,007.181 m for MODIS Sinusoidal ), X (Ae[-n,n]) the geo- 
graphic longitude in radian, 2 0 the longitude of the central 
meridian in radian, and (p {(p e [-n/2 t n/2]) the geographic lati- 
tude in radian (Seong et al., 2002). 

Mathematical transformation of UTM coordinates to coordi- 
nates in the MODIS Sinusoidal projection is achieved using PROJ.4. 
After converting the coordinates of four corners of a Landsat 
image into coordinates in the MODIS Sinusoidal projection, the 


B WRS-2 tiling system 


75° W 70° W 65° W 60° W 



A MODIS tiling system 



Fig. 3. Comparison of MODIS tiling system (A) and WSR-2 tiling system for Africa (B). (A) is made up of 36 paths and 18 rows. The paths numbered from 0 to 35, from 180E 
to 180W. The rows numbered from 0 to 17, from 90N to 90S. (B) for global is made up of 233 paths (numbered from 001 to 233, east to west, with path 001 crossing the 
equator at 64.60W) and 248 rows (numbered from 001 to 248, with row 001 starting at 80N and the numbering increases southward to a maximum latitude 81 S and then 
turns northward). 
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Fig. 4. Example of topological relationships between Landsat and MODIS tiles. A Landsat image can be located within a MODIS tile (A) or intersect with 2 (B), 3 (C), or 4 
(D) MODIS tiles. 


spatial coverage polygon of the Landsat image in the MODIS 
Sinusoidal projection can be reconstructed by connecting the four 
corners in a clockwise or anticlockwise order. When a Landsat 
image crosses the longitude line at 180E/180W, however, a spatial 
coverage polygon constructed this way will be wrong. To solve 
this problem, the area covered by the Landsat image is split into 
two parts along the 180E/180W longitudinal line. This way 
neither part crosses this longitude line, and the spatial coverage 
polygons of both parts in the MODIS Sinusoidal projection can be 
reconstructed following the procedure described earlier, which 
are eventually combined into a multipolygon geometry. 

After reprojecting the Landsat spatial coverage polygon using 
the MODIS Sinusoidal projection, the MODIS images needed to 
cover the target Landsat image can be identified based on their 
topologic relation to the Landsat image. Depending on the geo- 
graphic location, a Landsat image could be located within a 
MODIS tile, or interest with two to four Landsat tiles (see Fig. 4). 

4.3. Surface Reflectance Consistency Checking (SRCC) 

The Landsat-MODIS comparison process discussed in Section 3 
has been implemented in the SRCC module following OOP 
principles. OOP can facilitate the design and maintenance of 
model implementation especially when the system involves 
various data sources and subprocesses, or keeps growing by 
adding new algorithm components (Wolfsthal, 1994). In this case, 
OOP can benefit the system by organizing the algorithms in the 
comparison process and extending sample filtering and matrices 
algorithms in future. 

In OOP, the structure of object is defined by class. By inheriting 
a class (derived class) from another class (base class), the derived 
class reuses structure definition from the predefined base class 
(Thirunarayan, 1999). In addition to deriving attributes and 
functions from its base class, the derived class is also allowed to 
append more entities or even override the existing entities. 
Because the derived class contains everything that its base class 
has, it can act like its base class and be used to replace its base 
class without notifying the program units dependent on the base 
class. As such, inheritance can facilitate module extension by 
developing a derived class. Class could be designed only for 
extending purposes, and it is termed as abstract class opposed 
to concrete class. The abstract class has virtual functions which 
are required to be implemented by derived classes. 

The SRCC module consists of two class groups, i.e., framework 
classes and algorithm classes (see Fig. 5). The framework classes 
implement the key Landsat-MODIS comparison steps (see Section 
3.4) and basic elements associated with the comparison. Two 


abstract classes (i.e., SampleFiltering and IndictorCal class) are 
defined in the framework group to establish a way for adding 
algorithms classes into the comparison. The algorithm classes 
implement algorithms for sample filtering and metrics calcula- 
tion. Each algorithm class implements a particular algorithm and 
inherits one of the two abstract classes described earlier. 

4.3. 1 . Basic framework classes 

An object of Sample class represents one of the samples 
collected over the simultaneously acquired Landsat and MODIS 
images. In addition to the sample coordinate, the object includes a 
surface reflectance value pair and attributes related to the sample 
(e.g., MODIS zenith angle, MODIS cloud state, etc.). The surface 
reflectance value pair is collected from the MODIS and Landsat 
images at the sample coordinate, and the valid surface reflectance 
value range is from -100 to 16,000 (0.0001 scale factor) 
(Vermote and Kotchenova, 2008b). Samples with surface reflec- 
tance values out of the range are likely falling into the invalid 
image regions (e.g., image filling region, ETM+ scan line corrector 
(SLC) strips (Landsat 7 Science Data Users Handbook, 2009)), and 
are removed from the sample collection. The MODIS zenith 
value is retrieved from the 1 km resolution sensor zenith band 
(MODIS_Grid_l km_2D\SensorZenith_l). The MODIS cloud state 
is retrieved from the 1 km QA band (see Table 3). 

An object of the GeoRaster class encapsulates a geospatial 
image and associated operations. Its attributes include the geo- 
spatial coverage (extent), row and column number (size), and data 
blocks (tiles). Because loading large satellite images may consume 
lots of computer memory, data block is introduced to the system 
to allow loading a small portion rather than the whole image. 
Therefore, two functions are provided by the class to support 
reading specified pixel (readPixel) and data block (readTile). 

The Validation class implements the individual band compar- 
ison process (see Section 3.4) in its validationBand function. In 
addition, the validationlmage function performs comparison on 
the whole Landsat image by summarizing the comparison result 
for each band. 

4.3.2. Sample filtering classes 

The SampleFiltering class is an abstract class that defines the 
common attributes and functions needed for sample filtering 
algorithm classes. Each algorithm class is required to inherit the 
SampleFiltering class and implement its algorithm in the filter 
function. The algorithm class receives input samples from the 
samples attribute before calculation, and returns processed sam- 
ples after calculation. 
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Fig. 5. Classes defined for implementing the Landsat-MODIS comparison, including framework classes (top part) and algorithm classes (bottom part). Input and output 
parameters are wrapped as class attributes, and operations are wrapped as class functions. The attributes and functions of each class are divided into public or private 
groups according to their usage scope. Only the public scope attributes and functions are illustrated in the figure. 


Three algorithm classes were created: 

1. The CloudFiltering class detects cloud samples based on their 
MODIS cloud states. Samples with state other than “Clear” will 
be excluded from the comparison. 

2. The EquidirectionalFilter class excludes samples with MODIS 
zenith angle out of the Landsat zenith angle range ( ± 7.5°). 

3. The FIomogeneousFilter class detects and excludes the hetero- 
geneous observations of MODIS and Landsat following the 
rules discussed in Section 3.2.2. 


4.3.3. Metrics calculation classes 

The IndictorCal abstract class defines the common attributes 
and function for calculating the agreement metrics. Each algo- 
rithm class has to inherit the IndictorCal class and implement its 
algorithm in calculation function. The algorithm class receives 
collected samples through samples attribute before calculation, 
and returns agreement metrics through the indictor attributes 
after calculation. Two algorithm classes were created to imple- 
ment the linear regression ( LinearRegressionCal class) and RMSD 
calculation ( RMSDCal class), respectively. 

4.4. Outputs from Landsat-MODIS comparison 

Outputs from the Landsat-MODIS comparison include agree- 
ment metrics, scatter plots, and a shapefile showing the 


distribution of the agreement metrics for all Landsat images for 
a particular area. Agreement metrics include those discussed 
in Section 3.1 and are calculated for individual band (Landsat 
bands from 1 to 5, and 7) as well as by pooling all 6 bands 
together (see Table 4). Note that algorithm classes are extendable 
and more metrics can be calculated when new algorithms are 
added to the system. 

Scatter plots provide graphic views of the agreement metrics. 
A scatter plot is produced for each individual Landsat band, and 
an overall plot is generated by pooling all 6 spectral bands 
together (see Fig. 6). For continental to global scale applications 
that involve large quantities of Landsat images, a map showing 
the spatial distribution of the agreement metrics would allow 
quick identification of images having low agreement values. To 
meet such needs, the LMCCS also outputs the agreement metrics 
in shapefile format, which can be loaded and visualized as maps 
(see Fig. 9). 

4.5. System integration in an automated environment 

The LMCCS is developed using Java programming languages. 
Java has been widely supported by open communities. Imple- 
menting the system using Java can be benefited from millions of 
third part libraries and tools, especially the open geospatial 
libraries. Java is one of the OS independent programming 
languages, so the system can be used on any operation systems 
that support Java Virtual Machine (JVM), including UNIX/Linux 


M. Feng et al. / Computers & Geosciences 38 (2012) 9-22 


17 


Table 4 

Agreement metrics outputted for band and overall image from the Landsat-MODIS comparison. More metrics can be added when new metrics calculation algorithms are 
introduced into the system. 


Algorithms 

Agreement 

metrics 

Theoretical 

ranges 

Descriptions 

Linear regression 
calculation 

Slope (a) 
Offset (b) 
R 2 

[ — oo,oo ] 
[ — oo,oo ] 

[0,1] 

The slope and offset indicate the linear regression trend between Landsat and MODIS. For close trend, the slope is 
expected to be close to 1, and offset is expected to be close 0. 

The metric indicates the correlation between Landsat and MODIS. For good correlation, R 2 is expected to be close to 1. 

RMSD calculation 

RMSD 

[0, 00 ] 

The metric indicates the difference between Landsat and MODIS. The value is expected to be near 0 for little 
difference. 
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Fig. 6. Scatter plots for each individual Landsat band (6 small charts on right) and an overall plot for all 6 spectral bands (the big chart on left). MODIS and Landsat values 
are represented at x-axis and y-axis, respectively. The 1:1 line (red dash line) and the linear regression line (blue dash line) are added to each chart. The metrics are also 
labeled on the overall plot chart, including the linear regression (MODIS value, Landsat value, slope, and offset are labeled as x, y, a, and b respectively), R 2 , and RMSD. (For 
interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 


servers and personal computers. In order to facilitate automated 
use of the system, the system provides utility functions which can 
be executed as command lines. Thus, both users and other 
systems, such as LEDAP, will be able to use the system through 
those commands. Benefiting from the available Java-based Web 
server containers, e.g., Apache Tomcat (http://tomcat.apache.org/), 
GlassFish (https://glassfish.dev.java.net/), we believe that the system 
can be shared online with little modifications in future. 


5. Utilization of the LMCCS 

The overall goal of the LMCCS is to provide an efficient 
mechanism for verifying that surface reflectance products derived 
using Landsat images are consistent with those derived using 
MODIS data. High levels of consistency between the two products 
are indicators that both products are likely reliable and the 
reflectance values are comparable, and inconsistencies between 
Landsat and MODIS products as revealed by this system likely 
indicate quality problems with the Landsat products that need 


further investigation. Preliminary surface reflectance products 
derived from GLS 2000 epoch Landsat data have been used to 
evaluate the approach, and revealed problems in the products 
that were traced to incorrect version of the LEDAPS surface 
reflectance code and erroneous inputs. 


5.1. Incorrect version of the LEDAPS surface reflectance code 

An earlier version of the LEDAPS code obtained from NASA 
GSFC was used to produce a surface reflectance product using the 
GLS ETM+ image. The LMCCS revealed that the Landsat reflec- 
tance values were quite different from MODIS values in many 
places. Fig. 7 A shows the comparison results generated by the 
LMCCS for a cloud free image acquired over southern Turkey 
(WRS-2 path 174/row 34) on June 22, 2000. An investigation 
revealed that the module for retrieving calibration parameters 
from Landsat metadata files was designed for early versions of 
those files. After replacing that module with a new version 
designed to handle the metadata file format of the GLS data set, 
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Fig. 7. The LMCCS tool revealed significant disagreements between Landsat surface reflectance values calculated using an earlier version of the LEDAPS code and MODIS 
data (A). After an investigation (see text for details), a newer version of the LEDAPS code was obtained, which produced reflectance values highly consistent with MODIS 
data (B). In the example shown here, the Landsat 7 image was acquired on June 22, 2000, over southern Turkey (WRS-2 path 174/row 34). 
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Fig. 8. Low levels of agreements between Landsat and MODIS surface reflectance values (A) prompted an investigation into the inputs, which revealed that the input 
Landsat image (B) was corrupted. Another ETM+ image acquired at the same WRS-2 path/row but different date is shown in (C). Both (B) and (C) are shown with Landsat 
bands 4, 3, and 2 in red, green, and blue colors. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 


the LEDAPS code produced surface reflectance values highly 
consistent with the MODIS values (Fig. 7B). 

5.2. Erroneous GLS images 

The GLS data set consists of nearly 10,000 Landsat images 
produced using automated or semiautomated approaches, 
providing a near complete coverage of all land areas of the Earth 
(Tucker et al., 2004; Gutman et al., 2008). It is not clear whether 
each image has been checked to identify images that may have 
obvious quality problems. Since the LMCCS is designed for 
checking every Landsat image, it provides a mechanism for 
screening all GLS images to identify images that have obvious 
quality issues. Fig. 8A gives the Landsat-MODIS scatter plot 
generated by the LMCCS for a GLS 2000 acquired on May 12, 
2000, in Socotra, Yemen (WRS-2 path 159/row 51), which shows 
that there is essentially no correlation between Landsat and 
MODIS reflectance values. A visual check of the input image 


revealed that it was a corrupted image (Fig. 8B) compared 
to a normal image acquired at the same location (Fig. 8C). 
To our knowledge, this problematic image has never been 
identified after the GLS and its predecessor the GeoCover became 
available. 


5.3. MODIS surface reflectance calibration issue 

The LMCCS identified inconsistencies on band 7 of a Landsat 
surface reflectance image in the Sahara desert area. An investiga- 
tion revealed problematic values in its simultaneously acquired 
MODIS image. The issue has been reported to the MODIS team 
and confirmed to be a calibration error on band 7 of the MODIS 
image, which was acquired during the sensor’s early operation 
stage. Fig. 9 shows the Landsat image (Fig. 9A) acquired on April 
20, 2000, in Libya (WRS-2 path 181 /row 43) and its simulta- 
neously acquired MODIS image (Fig. 9B). 



M. Feng et al. / Computers & Geosciences 38 (2012) 9-22 


19 



Fig. 9. Low levels of agreements between Landsat and MODIS surface reflectance values at band 7 (A) prompted an investigation into the inputs, which revealed that the 
MODIS image band 7 has problematic calibrated pixels (illustrated in dark color) (B). (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this article.) 
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Fig. 10. R 2 distribution map of Landsat images for Southeast Africa. Most images have good correlation with the MODIS images. Using 0.8 as threshold, few images could 
be highlighted because of low R 2 . After manually checking the suspicious images, the image in Libya (A) has been confirmed to be a MODIS calibration error (see Fig. 9), the 
image on the northeast corner (B) has been confirmed to be a corrupted image (see Fig. 8), and the images near Madagascar Island (C) are caused by residual cloud. 


5.4. Consistency check from regional to global scales 

To improve the efficiency of checking large quantities of 
Landsat images needed to cover a large area, the LMCCS also 
outputs the agreement measures into a shapefile, which can be 
used to create maps showing the spatial distribution of those 
agreement measures. By visually inspecting these maps, areas 
where Landsat and MODIS reflectance values are inconsistent can 


be identified quickly. Fig. 10 shows a map of the R 2 value 
generated from the Landsat-MODIS scatter plots for 723 GLS 
images covering most areas in eastern Africa. It shows that of the 
723 images evaluated here, only four of them had R 2 values less 
than 0.8. Further investigations on the four images revealed that 
one of them had a problem with the original input image (see 
Section 5.2), another one caused by problematic calibration of 
MODIS (see Section 5.3), and the other two had very small 
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portions of land area. These two images have lots of clouds over 
water, many of which may not have been flagged properly and 
contributed to the low correlations between the Landsat and 
MODIS data. 


6. Discussion and conclusions 

Advances in developing automated atmospheric correction 
algorithms coupled with available computing and storage capa- 
cities have made it feasible to produce large quantities of surface 
reflectance products using Landsat images (Masek et al., 2006). 
During the process of converting raw satellite radiometry to 
surface reflectance, however, errors could be introduced from 
many potential sources, including incorrect input parameters, 
corrupted input images, bugs in the atmospheric correction 
algorithm or use of an incorrect version of the algorithm, and 
operator errors (Roy et al., 2002). To rule out the impact of such 
potential error sources, it is necessary to have a sanity check on 
the surface reflectance values derived from each Landsat image. 

Because MODIS surface reflectance products have been vali- 
dated comprehensively and are available globally since 2000, they 
can be used as a reference to validate the radiometric quality of 
surface reflectance products derived from essentially all Landsat 
images acquired during the MODIS era. Such a quality assessment 
approach is especially valid for Landsat 7 images because the 
MODIS on the Terra platform and the Landsat 7 share the same 
orbit and are only half an hour apart, and each of the 6 Landsat 
spectral bands overlaps with a MODIS band. In general, good 
agreements between MODIS and Landsat surface reflectance 
values can be considered indicators of the reliability of the 
Landsat products, while disagreements may suggest potential 
quality problems with the surface reflectance products that need 
to be further investigated. However, a number of factors can lead 
to low agreement between Landsat and MODIS reflectance values 
even when both products are perfectly valid, including pixel 
misalignment due to improper handling of pixel and image 
geometry, fine scale spectral variations coupled with uneven 
weightings of sensor’s point spread function, and significant 
changes in atmospheric conditions that may occur between the 
overpass of the Landsat 7 and the Terra MODIS. 

The LMCCS is designed to perform automatic comparison of 
Landsat and MODIS surface reflectance data while eliminating or 


minimizing the impact of the factors described above on the 
derived agreement metrics. In this tool, MODIS and Landsat pixels 
are aligned properly using a widely used reprojection tool and by 
tracking the geometry of pixels in both data types accurately. The 
impact of fine scale spectral variations and sensor’s point spread 
function is minimized by only using samples from homogeneous 
regions in calculating agreement metrics. For Landsat images 
having minimal cloud cover, which in theory should be case for 
the GLS images, the implemented cloud screening mechanism 
that only considers clouds in the MODIS data can effectively 
reduce the impact of changes in cloud location and/or radiometry 
that may have occurred between the overpasses of the Landsat 
7 and the Terra MODIS (Fig. 11). For more general use of the 
LMCCS regardless of cloud cover in the Landsat images, the cloud 
screening mechanism should consider clouds in both the Landsat 
and the MODIS data. 

To improve the robustness of the LMCCS in an automated 
environment, a number of issues have been considered in devel- 
oping the LMCCS. The geospatial MODIS Image Metadata Data- 
base allows quick identification of simultaneously acquired 
MODIS images needed to cover any Landsat image. MODIS data 
gaps resulting from various practical reasons have been identified 
properly, and the database can be updated as more MODIS data 
are made available. Many reprojection tools have difficulties in 
handling images that encompass the 180E/180W longitudinal 
line. This situation has been handled properly in the LMCCS. 

To maximize its portability, the LMCCS is developed using Java 
and open-source libraries. It does not require any commercial 
software and can be deployed on any platform that supports JVM. 
In 2008, NASA and USGS implemented a free Landsat Data 
Distribution Policy that provides level 1 terrain corrected data 
for the entire U.S. Landsat archive, and the needs of quality 
assessment are increasing while more Landsat surface reflectance 
are expected to be produced from the free Landsat data. The 
Java-based LMCCS has great potential to be modified to be an 
online system, providing operational quality assessment ability to 
more users through the Internet. OOP principles are followed in 
order to make this tool more flexible for future expansion, 
including use of new filters and sampling schemes, adding new 
agreement metrics, and adaption for use with other Landsat class 
instruments such as the China-Brazil Earth Resources Satellite 
(CBERS), the Linear Imaging and Self-Scanning sensor (LISS-III) 
and Advanced Wide Field Sensor (AWiFS) onboard the Indian 
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Fig. 11 . Clouds in an image (C) can significantly reduce the agreement values derived from a Landsat-MODIS comparison (A) as they can move a lot between the 
overpasses of the Landsat 7 and the Terra MODIS. If the Landsat reflectance values are reliable, the agreement values are improved greatly after the clouds are excluded 
using the cloud screening mechanism implemented in the LMCCS (B). The image in (C) is shown with Landsat bands 4, 3, and 2 in red, green, and blue colors. Most clouds 
appear bright white in the lower left corner of this image. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this 
article.) 
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Remote Sensing Satellites (IRS), and the Advanced Spaceborne 
Thermal Emission and Reflection Radiometer (ASTER) onboard 
NASA’s Earth Observing System (EOS) Terra spacecraft. This 
approach also has the potential to be adapted for quality assess- 
ment of surface reflectance produced from next generation Land- 
sat sensors, the LDCM. Due to orbital differences, for some of 
these instruments, including the Landsat 5, the Terra MODIS may 
not provide a daily image that is acquired on the same day as the 
image to be assessed and has illumination and viewing geometry 
similar to that image (Landsat 7 Science Data Users Handbook, 
2009). In such cases, the MODIS Nadir BRDF-Adjusted Reflectance 
(Schaaf et al., 2002) product may be the best alternative that can 
be used in the place of the MOD09GA data. 
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